Viscoelastic subdiffusion: from anomalous to normal 
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We study viscoelastic subdiffusion in bistable and periodic potentials within the Generalized 
Langevin Equation approach. Our results justify the (ultra)slow fluctuating rate view of the corre- 
sponding bistable non-Markovian dynamics which displays bursting and anti-correlation of the res- 
idence times in two potential wells. The transition kinetics is asymptotically stretched-exponential 
when the potential barrier Vb several times exceeds thermal energy ksT (Vo ~ 2 -j- 10 ksT) and 
it cannot be described by the non-Markovian rate theory (NMRT). The well-known NMRT result 
approximates, however, ever better with the increasing barrier height, the most probable logarithm 
of the residence times. Moreover, the rate description is gradually restored when the barrier height 
exceeds a fuzzy borderline which depends on the power law exponent of free subdiffusion a. Such 
a potential-free subdiffusion is ergodic. Surprisingly, in periodic potentials it is not sensitive to the 
barrier height in the long time asymptotic limit. However, the transient to this asymptotic regime 
is extremally slow and it does profoundly depend on the barrier height. The time-scale of such 
subdiffusion can exceed the mean residence time in a potential well, or in a finite spatial domain by 
many orders of magnitude. All these features are in sharp contrast with an alternative subdiffusion 
mechanism involving jumps among traps with the divergent mean residence time in these traps. 

PACS numbers: 05.40.-a, 82.20.Uv, 82.20.Wt, 87.10.Mn, 87.15. Vv 



I. INTRODUCTION 

Multifaceted anomalous diffusion attracts ever increas- 
ing attention, especially in the context of biological appli- 
cations. For example, diffusion of mRNAs and ribosomes 
in the cytoplasm of living cells is anomalously slow [l[, 
large proteins behave similarly [2|, y, \% |5| . Even intrinsic 
conformational dynamics of the protein macromolecules 
can be subdiffusive @, 0, H, i, M, ED, El El, EI El- 
There is a bunch of different physical mechanisms and 
the corresponding theories attempting to explain the ob- 
served behaviors, from spatial and/or time fractals, in- 
fluence of disorder, cluster percolation, etc., to viscoelas- 
ticity of coinplex media [ij, E3, El, 111 HE HI, H [H, 
HI, HH, HE E3, Hl| ■ I n particular, molecular crowding can 
be responsible for the viscoelasticity of dense suspensions 
like cytosol of bacterial cells lacking a static cytoskeleton 
@, US Ell- The state of the art remains rather perplexed, 
offering cardinally different views on the underlying phys- 
ical mechanisms, as we clarify further with this work. 

One physical picture reflects a set of the traps (possi- 
bly dynamical (27() where the diffusing particle stays for 
a random time r». The mean residence time (MRT) in 
traps should diverge [H, Ell for the diffusion to become 
anomalously slow, i.e. with the position variance growing 
sublinearly, (5x 2 (t)} ~ t a , with < a < 1. This stochas- 
tic time-fractal picture became one of the paradigms in 
the field [l7[- It can be also related to averaging over 
static, or quenched disorder [T3, H(|. Such a continu- 
ous time random walk (CTRW) among traps has infinite 
memory even if the residence times in the different traps 
are not correlated. The memory comes from the non- 
exponential residence time distributions in traps. 

A quite different physical view was introduced by Man- 
delbrot et al. 13111 with the fractional Brownian mo- 



tion (fBm). Here, the standard Gaussian Wiener process 
with independent increments is generalized to incorpo- 
rate the statistical dependence of increments. The Gaus- 
sian nature remains untouched, but the increments can 
be either positively, or negatively correlated over an infi- 
nite range. Positive correlations (persistence) lead to su- 
perdiffusion. If correlations are anti-persistent, i.e. given 
a positive increment, the next one will, with greater prob- 
ability, be a negative increment and vice versa, then a 
subdiffusive behavior can result. This idea does not im- 
ply that the residence time in a finite spatial domain 
diverges on average. Here roots the cardinal difference, 
in spite of some superficial similarities, between the fBm- 
based and the CTRW-based approaches to subdiffusion. 

FBm emerges naturally, e.g., in viscoelastic media as 
one of the best justified models. Indeed, let us start 
from a phcnomcnological description of viscoelastic forces 
acting on a particle moving with velocity x(t) in some 
time-window [0,t): 

F v - el (t) = - f r)(t-t')x(t')dt'. (1) 
Jo 

Clearly, for a memory-less linear frictional kernel with 
rj(t) = 2rj8(t) on the particle acts a purely viscous Stokes 
friction force, F v = —r\x. If memory does not de- 
cay, 77(f) = r\ = const, then the force is quasi-elastic, 
F e i(t) = —rj[x{t) — x(0)] (cage force). A popular model 
of viscoelasticity introduced by Gemant (32|, which in- 
terpolates between these two extremes, corresponds to 
rj{t) = n a t~ a /T(l - a) with < a < 1 (T(x) is 
the gamma- function). Remarkably, this model yields 
the Cole-Cole dielectric response for particles trapped in 
parabolic potentials [33ll34|. which is frequently observed 
in complex media. For a small Brownian particle of mass 
to, one must take into account unbiased random forces 



■2 



acting from the environment (Langevin approach). 
Then, the linear friction approximation combined with 
the symmetry of detailed balance fixes the statistics of 
the stationary thermal random forces to be Gaussian 
[35l ]. Moreover, the fluctuation-dissipation theorem dic- 
tates that the stationary autocorrelation function of the 
random force, temperature T, and the memory kernel are 
related 1361] : 



(mm)=k B Tr 1 (\t~t'\) 



(2) 



For Gemant model, is the fractional Gaussian noise 
(fGn) [TH, [31], [3ij]. Altogether, the motion in potential 
V(x, t) is described by the generalized Langevin equation 
(GLE) 



mx + n(t - t')x(t')dt' + 



dV(x,t) 
dx 



m . (3) 



Importantly, this GLE can also be derived from the me- 
chanical equations of motion for a particle interacting 
with a thermal bath of harmonic oscillators, i.e. from 
first principles. This statistical-mechanical derivation 
[H, [37], [H, H^] involves the spectral density J(u>) of 
bath oscillators. It is related to the spectral density 
of thermal force, S(oj) = 2 L (£(t)£(0)) cos(ut)dt, as 
S(w) = 2k B TJ(ui)/uj. With J(uj) = r] a ui a , the so-called 
Ohmic case of a = 1 corresponds to viscous Stokes fric- 
tion and normal diffusion. The sub-Ohmic, or fracton 
thermal bath [H, [H with < a < 1 and noise 
spectrum of random force corresponds to the above Ge- 
mant model of viscoelasticity. It yields subdiffusion in 
the potential- free case [1^, [4(| . The velocity autocorre- 
lation function is then negative (except of origin), being 
the reason for the anti-persistent motion. The physical 
origin of this feature is that the elastic component of the 
viscoelastic force opposes the motion and ever tries to 
restore the current particle's position. Moreover, in the 
inertialess limit (m — > 0) the solution of GLE is fBm 
with the coordinate variance (Sx 2 (t)) — 2K a t a /L(l + a) 
[H, [4l[ and the subdiffusion coefficient K a obeying the 
generalized Einstein relation K a = kBT/r/ a [42, 43]. This 
way, anti-persistent subdiffusive fBm emerges from first 
principles within a physically well grounded, but approxi- 
mate description. It corresponds also exactly to the diffu- 
sion equation with a time-dependent diffusion coefficient 
D(t) = K a /[T(1 + a)i 1_Q; ] which is frequently used 
to fit experiments in viscoelastic and crowded environ- 
ments, see e.g. in pi. 15, |29| . 



II. THEORY 

Anomalous escape (rate) processes and spatial sub- 
diffusion in periodic potentials represent within the 
GLE description a highly nontrivial, longstanding chal- 
lenge. Even the corresponding non-Markovian Fokker- 
Planck description is generally not available, except 
for the strictly linear and parabolic potentials [441 ]. 
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FIG. 1: (Color online) Frictional power law memory kernel 
(in units of r) a ) and its two different approximations versus 
time t (arbitrary units) for a — 0.5 and vq = 10 3 . Notice that 
the approximation with b = 2 practically coincides with the 
exact kernel in this plot and the choice b = 10 is also a very 
good one in spite of logarithmic oscillations which are barely 
seen. Inset magnifies a part of plot. 



To get insight into the physics of such processes, it 
is convenient to approximate 1// 1_Q noise by a sum 
of independent Ornstein-Uhlenbeck noise components, 
£(t) = X^o 1 CiW> w ith t ne autocorrelation functions, 
(Ci(t)Cj(t')) — ksTrjiSij exp(— j/j|t — The correspond- 
ing memory kernel is accordingly approximated by a sum 
of exponentials, 



JV-l 

77 (t) = 22, ??iexp(-t^), 

i=0 



(4) 



where Vi = vo/b z is the inverse autocorrelation time of 
the i-th component and r)i = (r] a /T(l — a)) C a {b)vQ jb ta is 
its weight. Furthermore, vq presents the high-frequency 
cutoff of £(i), b is a dilation (scaling) parameter, and 
C a (b) is a numerical constant. The low- frequency noise 
cutoff is uj c = vo/b 1 * -1 . It is worth to mention that such 
cutoffs emerge for any 1// noise on the physical grounds 
l45ll . For a — 0.5, which is of experimental interest [H, 
Il4l ]. the choice of b = 2 (i.e. octave scaling) and N — 
64 with vq = 10 3 (in arbitrary units) and 61/2(2) = 
0.389 allows one to fit perfectly the power law kernel 
in the range from t = 10~ 3 till t — 10 15 , i.e. over 18 
time decades. The choice of b = 10 (i.e. decade scaling) 
with Ci/2(10) = 1.3 provides also an excellent fit over 
15 time decades from t = 10~ 3 till t = 10 12 with N = 
16. The numerical advantage of larger b is that one can 
use smaller Markovian embedding dimension D — N + 
2. These two approximations to the exact power law 
memory are shown in Fig. [TJ The approximation with 
b = 10 displays logarithmic oscillations [l7| which are 
barely seen in this plot and make a little influence on the 
stochastic dynamics, see in Fig. [5] 

Free subdiffusion holds until the time scale of l/u) c 
which can be very large. The idea of such a represen- 
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tation of a power law dependence is rather old [17J, [46J , 
being also habitual in the 1 // noise theory [45| . The cor- 
responding power spectrum is approximated by a sum 
of Lorentzians, S(uj) — 2k B T ^\ Jl+l't ■ Every station- 
ary noise component is asymptotic (t — > oo) solution 
of Q(t) = -ViQ(t) + y/27]ii/ik B T£i(t), where are 
independent white Gaussian noises with unit intensity, 
(£,i(t)£,j(t')) = SijS(t—t'). Furthermore, the particle must 
act back on the source of noise in order to have the FDT 
relation © satisfied. This yields the following D = N+2 
dimensional Markovian embedding of the non-Markovian 
GLE stochastic dynamics in Eq. ([3]) with kernel ([!]): 

x = v , (5) 
mv = -^p^+J^uM, (6) 

X i=0 

■iii = -rjiV - v^i + \j2vir\ik B Tii(t) . (7) 

Initial itj(O) have to be sampled independently from un- 
biased Gaussian distributions with the standard devia- 
tions <7j = y/ksTrji [47]. Under this condition, it is easy 
to show that Eqs. ©-© are equivalent to the GLE © 
with kernel © under FDT relation ©. Notice that itj(t) 
are auxiliary mathematical variables and should not be 
interpreted as (scaled) coordinates of some physical par- 
ticles. The embedding ©-([7]) can be also derived from a 
more general scheme in Ref. [47i j . 

The numerical simulations of Markovian dynamics in 
Eqs. ©-([7]) below were done with stochastic Euler and 
Heun algorithms [48| using different random number gen- 
erators. The results are robust. The simulations have 
been checked against the exact analytical results avail- 
able for the potential-free case and parabolic potentials. 
The both above embeddings with b = 2, N = 16 and 
b = 10, N = 64 yield practically the same results within 
statistical errors. However, the simulations with N = 16 
requires much less computational time. Therefore, we 
preferred the latter D = 18-dimensional embedding in 
most simulations following a "rule of thumb" [49( : a neg- 
ative power law extending over n time-decades can be 
approximated by a sum of about n exponentials. The 
quality of this approximation along with numerical errors 
is discussed in the Appendix A and Fig. [5] by making 
comparison of the numerical Monte Carlo results with 
the numerically exact solution of the free subdiffusion 
problem. The numerical error is mostly less than 3% for 
free subdiffusion in this work, whereas the theoretical er- 
ror incurred by the 16-exponential approximation of the 
power law memory kernel is mostly less than 1%, cf. Fig. 
[2j Clearly, it makes no sense to approximate the kernel 
better, if no more than n = 10 4 trajectories are used in 
the ensemble averaging. 

The chosen Markovian embedding of non-Markovian 
GLE dynamics with r){t) in Eq. © is mathematically, 
of course, not unique. Another embedding was proposed 
in Ref. [501 ] and infinitely many different embeddings of 
one and the same non-Markovian dynamics are in fact 
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FIG. 2: (Color online) (a) Stochastic simulations of Eqs. 
©-0 with embedding dimension D = 18 vs. exact solu- 
tion of GLE <(3j for free subdiffusion. Coordinate is scaled 
in some arbitrary units L and time is scaled in units of 
T0 = (L 2 77 a /fcflT) 1 ' a ; a = 0.5 and the damping parameter 
rp = r^/Tbai = 10 with r ba i = L/v T and v T = ^Jk B T/m. 
Particles are initially localized at x = with Maxwellian dis- 
tributed velocities. The number of particles n = 10 4 is used in 
simulations (stochastic Heun algorithm, time step At = 10~ 4 ) 
for the ensemble averaging. The time averaging for a single 
trajectory is done using Eq. ()A5|I . The agreement proves er- 
godicity. Inset in compares two numerically exact solutions 
of GLE, one with the strict power law kernel and one with 
its 16-exponential approximation, (b) Relative numerical er- 
ror of our simulations (versus the exact solution) is presented. 
The relative deviation of the solution with the approximate 
memory kernel from the solution with exact power law kernel 
is also presented. Notice the occurrence of logarithmic oscil- 
lations within a less than 1% error margin, which, however, 
do not play any essential role. 



possible [47|. However, our simplest scheme allows to 
contemplate straightforwardly the view of anomalous es- 
cape processes as rate processes with dynamical disorder 



4 



A. Non-Markovian rate theory and beyond 

We consider now stochastic transitions in a paradig- 
matic bistable quartic potential V(x) = V (l — x 2 /xq) 2 
with minima located at x m - m = ±x and the barrier 
height Vq. The question is: Which is the statistical dis- 
tribution of the residence times in two potential wells and 
the escape kinetics? This is a long-standing problem for 
a general memory friction. Since the effective friction is 
sufficiently strong (the memory friction integral diverges) 
one can tentatively use a prominent non-Markovian rate 
theory (NMRT) result [HI [H, which is a general- 
ization of the celebrated Kramers rate expression [57| . It 
assumes asymptotically an exponential kinetics for the 
survival probability in one well, P(t) ~ exp[— R(p)t], 
with the non-Markovian rate 

RQM) = Kfa)^exp(-l3V ). (8) 

In Eq. ©, uj = y/V" (x min )/m = yj8V /(mxl) is the 
bottom attempt frequency, cxp(— /3Vq) is the Arrhcnius 
factor, /3 = l/(fcgT) is the inverse temperature, and 

= - < 1 (9) 

is the transmission coefficient. It invokes the effective 
barrier frequency /i given by the positive solution of equa- 
tion 

M 2 +m{n)/m = w 2 , (10) 

where f](s) — exp(— st)rj(t)dt is the Laplace- 

transformed memory kernel, and uib = \J\V" (Q)\/m = 
ujq I is the (imaginary) barrier top frequency in the 
absence of friction. We focus below on the case of suffi- 
ciently high barriers, where the Arrhenius factor is small, 
exp(— (3Vo) <C 1. Clearly, a good single-exponential 
kinetics with exponentially distributed residence times, 
tp(r) = —dP(r)/dr, can only be valid for such potential 
barriers, even in the strictly Markovian case. 

However, how high is high? Could asymptotically ex- 
ponential kinetics be attained for the viscoelastic model 
considered at all? Very important is that the relaxation 
within the potential well is ultraslow and this fact seems 
to invalidate the non-Markovian rate description gener- 
ally [4l|]. To understand this, let us neglect formally for 
a while the inertia effects, m — ► 0. Then, the strict power 
law kernel corresponds (in parabolic approximation) to 
the relaxation law (Sx(t)) = 5x{Q)E a [—{t/T r ) a ] with 
the anomalous relaxation constant r r = [ri a x 2 t /(8Vo)] 1 ^ a , 
where E a (x) = yi ^°_n JC n /r(l + an) is the Mittag-Leffler 
function [IH, l4l| . Asymptotically, (Sx(t)) oc (t/r r )~ a , 
being initially a stretched exponential. Precisely the 
same relaxation law holds also for the CTRW subdif- 
fusion in the parabolic well (43| which (along with other 
similarities for the potential-free case) gave grounds to 
believe that these two subdiffusion scenarios are some- 
how related, or similar. From the fact of ultraslow re- 
laxation, it is quite clear that there cannot be a rate 



description even for appreciably high potential barriers, 
until the relaxation time within a potential well becomes 
negligible as compared with a characteristic time of es- 
cape. It worth to mention here that the non-Markovian 
rate theory approach yields a finite rate always, even 
for the strict power-law kernel [58[, Rq = R(u!r ) = 
^ &) exp(-/?y )/(\/27r), where J r b) = [4V / '{n a xl)) l / a is 
solution of Eq. (fTOf) for m — > 0. 

This cannot be, however, always correct. Indeed, let us 
introduce ad hoc a variable small-frequency cutoff v c such 
that fj(s) becomes fj(s) = n a (s + v c ) a ■ Then, choosing 
self- consistently v c — R{p) in Eqs. (|5|)- (fTU]) (for m — ► 0) 
one can show that the corresponding fi becomes modified 
as fi = U b) [l + exp(-/3Vb)/(\/27r)] 1 / Q - 1 . From 

this we conclude that the non-Markovian rate expression 
i?(/i) is practically not affected by such a cutoff when 
exp(— (3Vq) <C 1. This does not mean, however, that 
all the slowly fluctuating noise contributions with z/j < 
Rq can be simply neglected. They lead, in fact, to the 
fluctuating rate description invalidating thereby the non- 
Markovian rate picture. 

B. Fluctuating rates: simplest approximation 

The idea is to divide all the noise components Q into 
the two groups, = + £, s (t)- the fast noise 

= Tli<i which contributes to the "frozen" 

non-Markovian rate R(fJ-), and the slow modes which 
constitute the slowly fluctuating random force £ s (i) = 
Tlii>i d{t)- The separation frequency z/; c is chosen such 
that Vi < Vi a < R(n). It depends on the ratio of bar- 
rier height and temperature, as well as a. Furthermore, 
let's assume that for the slow Ui-modes in Eq. ([7]) one 
can approximately replace v(t) by its average zero- value. 
This is a reasonable approximation because of the dy- 
namics of v(t) is fast on that time scale. Then, the cor- 
responding equations for Ui(t) decouple from the parti- 
cle dynamics, Uj — ► and the corresponding stochas- 
tic modes can be considered just as an external random 
force. The fast noise agitates the particle trapped (other- 
wise) in the potential wells leading to the escape events. 
To a first approximation, one can regard the slow noise 
£ s (£) be quasi- frozen on the time scale of such escape 
events. Then, for high barriers /3Vq 3> 1 one can use a 
two-state approximation for the overall kinetics with the 
non-Markovian rate R(fj,*,^ s ) slowly driven in time by 
£, s (t) . This slow stochastic force is in fact also power-law 
correlated. Thus, we are dealing with a typical problem 
of non-Markovian dynamical disorder . Some insight 
can be obtained by using the quasi-static disorder ap- 
proximation [13, [U, [E^, [6(| for the averaged kinetics, 

w(^)exp[-R h2 (^,^)t}, (11) 

-00 

where JZi 2(^*1 Ca) are the non-Markovian rates 
for a quasi-static biasing force £ s distributed 



with the Gaussian probability density w(£ s ) = 
l/(^/2na s ) exp(— £ 2 /(2cr 2 )) and variance a 2 = k B Tr) s , 
where r] s = J2i>i Vi- The calculation of a 8 for a = 0.5 
(or larger) shows that the bias fluctuations are suffi- 
ciently small for (3V$ > 2, so that the approximation 
#i,2(m*>£s) ~ R(p*) exp[±^ s x /(fcsT)] can be used. 
Here, we just assume that the rms of potential barrier 
modulations ±a s xo is small against Vq. Since the influ- 
ence of slow modes on the effective barrier frequency /i 
is exponentially small for high barriers (see above), one 
can replace R(n*) with R(n). This finally yields 

POO 

Pi, a (t)« / exp[-R(n)ty]W{v)dy, (12) 
Jo 

where W(y) = l/(V2wdy) exp[- ln 2 (y)/(2d 2 )] is the 
probability density of log-normal distribution with width 
d = \J i] s Xq/ (ksT). The corresponding mean residence 
time (MKT) and the relative standard deviation, 5a = 
V(r 2 )-{r) 2 /{r), are: 

<r) = J R- 1 (/i)exp(d 2 /2), (13) 

5a = V2cxp(d 2 )-1 = y/2[(T)R(/i)] 2 - 1. (14) 

To characterize non-Markovian kinetics, one can intro- 
duce also a time-dependent rate k(t) = —d\nP{t)/dt 
which decays asymptotically to zero for any finite width 
d within this approximation. 

The resulting physical picture becomes clear: Fast 
Ornstein-Uhlenbeck components with > i?(/z) partici- 
pate in forming the non-Markovian rate -R(/-t), while the 
slow ones lead to a stochastic modulation of this rate 
in time. This implies the following main features which 
are confirmed further by a numerical study: (i) Both 
the mean residence time in a potential well, and all the 
higher moments exist, (ii) Anti-correlations between the 
alternating residence time intervals in the potential wells 
emerge along with a profoundly bursting character of the 
trajectory recordings, cf. Fig. [3J Indeed, during the time 
of a quasi-frozen stochastic tilt many transitions occur 
between the potential wells. The shorter time in the 
(temporally) upper well is followed by a longer time in 
the lower well. This yields anti-correlations (cf. Fig|4}. 
Moreover, many short-living transitions into the tempo- 
rally upper well occur, which appears as bursting (cf. Fig. 

The subsequent sojourns in one potential well are also 
positively correlated for many transitions (not shown), 
(iii) The escape kinetics is clearly non-exponential, see 
in Fig. [5] and below, (iv) The corresponding power 
spectrum of bistable fluctuations has a complex struc- 
ture with several different 1// noise low-frequency do- 
mains (cf. Fig. [6]); (v) The higher is potential barrier, 
the smaller is the rms of slow rate fluctuations. The 
last circumstance implies that for very high barriers the 
exponential escape kinetics with non-Markovian rate in 
(|8l)-lflQ|) will be restored, cf. Fig. [H For small a, this 
however can require very high barriers and be practically 
unreachable. 
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FIG. 3: (Color online) A sample trajectory of bistable 
transitions for (3Vo = 6. Time is given in units of r^ 6 ' = 
D7a£o/(4Vo)] 1//a and coordinate in units of xq. The residence 
time distributions are extracted by using the two thresholds 
(red lines). Stochastic Heun algorithm with the time step 
At = 10 -3 and Mersenne Twister pseudo-random number 
generator, combined with the Box-Muller algorithm are used. 
Markovian embedding with N — 16 is used as described in 
the text, a = 0.5, r = 10, b = 10, u = 1000. 



C. Numerical results 

Let us compare now these theoretical predictions with 
numerical results. The time is scaled in this section in 
the units of rl h) = 1/W ] = [rj a xl/(AV )] 1/a , which is the 
anomalous relaxation constant for the inverted parabolic 
barrier in the overdamped limit, and the role of inertial 
effects is characterized by the dimensionless parameter 
r = uibTr ■ The used r — 10 corresponds to the over- 
damped limit in the case of normal diffusion. For the 
used 16-exponential approximation of the memory ker- 
nel, it yields ijl « 0.999 in Eq. (fTUj) which is very close 
to /i = 1 corresponding to the formal overdamped limit, 
to — > 0, with the transmission coefficient n = l/r = 0.1. 
Despite this fact, some inertia effects for the intra- well 
relaxation dynamics are still present. Generally, it is im- 
portant to include such effects for a power law memory 
kernel (olT |. We performed simulations of very long tra- 
jectories (from 5 • 10 3 to 10 6 transitions between wells) 
achieving statistically trustful results in each presented 
case. A sample of stochastic trajectory for (3Vo = 6 is 
shown in Fig. [3] The bursting character is clear [5!|, 
indicating also slow tilt fluctuations. 

To extract the residence time distributions (RTDs) 
in the wells, Vi(ti) and 1^2(^2), and their joint distri- 
bution ■0(ti,T2), two thresholds were set at the min- 
ima of the potential wells, cf. Fig. [3] Fig. Q] displays 
the the joint distribution ^(ln(ri), ln(ra)) of the loga- 
rithmically transformed residence times for f3Vo = 2. 
Two facts are self-evident: (1) the transformed distribu- 
tion is not sharply peaked and spreads over several time 
decades; (2) the subsequent residence times in two poten- 
tial wells are significantly anti-correlated. The normal- 
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FIG. 4: (Color online) Two-dimensional distribution of the 
log-transformed residence times in the potential wells for 
PVo = 2. The green symbol corresponds to the non-Markov 
rate theory result (for the studied 16-exponential expansion 
of the power-law memory kernel). The orientation and the 
structure of the plateau of maximal probabilities manifest 
anti-correlation of the residence times in two potential wells. 
This feature is completely beyond the non-Markovian rate 
theory. Statistical software R 62] is used for data analysis 
and to produce the plot. 



ized covariance between n and Ti is c(ti,T2) ~ —0.116, 
and between the logarithmically transformed variables 
c(lnTi,lnT2) k, —0.19. The mean residence time is ap- 
proximately (ti^) ~ 52 with the relative standard devia- 
tion (5(71,2 = \j (T~i t2 ) ~ ( t i,2) 2 /(ti,2) ~ 1-69, whereas the 

non-Markovian rate theory yields tnm = 1 / JS(/Lt) ~ 32.9. 
This value essentially underestimates MRT, but it lies 
not far away from the extended region of most probable 
In Ti. 2 (see the green symbol in Fig. |4|). Furthermore, 
the distribution of the residence times in each potential 
well, tp( T ) — ~P( T )i is profoundly non-exponential, with 
a complex kinetics being mostly stretched-exponential, 
P(t) ~ exp[— (r/r 7 ) 7 ]. The power 7 slightly varies in 
time and reaches asymptotically 7 « 0.65, as indicated 
by a straight line trend for — ln[P(r)] on the doubly- 
logarithmic plot in Fig. [5] Generally, the asymptotic 
value of 7 is bounded as a < 7 < 1 and depends on 
the ratio Vb/fcsP- The formally defined time-dependent 
non-Markovian rate decays to zero as k(t) oc 1/t 1 1 
and the corresponding power spectrum of fluctuations in 
Fig. [6] displays a complex 1 / / 7 -noise pattern with the 
same 7 ~ 0.65 at lowest frequencies. Overally, the non- 
Markovian rate theory approach is clearly not applicable 
for such a barrier. 

The qualitatively similar features remain also for some 
higher potential barriers (or lower temperatures), e.g. for 
(3Vq = 6. In this case, numerically (r x 2 ) ~ 2710, and 
<5<7i.2 ~ 1.41, whereas the non-Markovian rate theory 
yields r NM « 1794 with ln(r NM ) « 7.49. This NMRT re- 



FIG. 5: (Color online) Survival probability in one potential 
well versus time (in units of Tr = [n a x^/{4Vo)] 1 ^ a ) for /3Vb = 
2. The asymptote is stretched-exponential with 7 w 0.65. Fit 
with Eq. is done using d = 0.81 and R = 2.67 • 10" 2 . 

These values were derived from the numerical data using the 
first two moments of the numerical distribution and Eqs. (|13|l , 
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FIG. 6: (Color online) The low-frequency part of the power 
spectrum of fluctuations for /3Vo = 2. The 1/f feature 
for smallest frequencies is defined chiefly by the stretched- 
exponential asymptotics of the survival probability, cf. Fig. [3] 
The correlogram method along with the SSA-MTM Toolkit 
for Spectral Analysis [63j is used to produce the spectrum. 



suit compares, however, now well against the most prob- 
able ln(ri 2) m Fig- This provides one of important 
results: Even if the non-Markovian rate theory is still not 
applicable, it can predict remarkably well the most prob- 
able logarithm of residence times. The kinetics remains 
asymptotically stretched-exponential even for such a high 
barrier with 7 increased to 7 w 0.76 (cf. Fig. [8|). How- 
ever, the region of most probable lnri^ shrinks further 
with increasing (3Vq and the non-Markovian rate theory 
describes ever better both the most probable lnri ,2, and 
the (logarithm of) mean residence time which start to 
merge as it should be for a single-exponential RTD. Al- 
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FIG. 7: (Color online) Two-dimensional distribution of the 
log-transformed residence times in the potential wells for 
PVo = 6. The non-Markovian rate theory result (green sym- 
bol) yields the most probable value of ln(ri,2). Kinetics is, 
however, clearly non-exponential, even asymptotically. 
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FIG. 8: (Color online) Survival probability in one potential 
well versus time (in units of t^) for /3Vo = 6. Asymptotics 
is stretched exponential with 7 ~ 0.76. Initial kinetics is 
stretched exponential with 7 m a. Fit with Eq. (|12[1 is done 
using d = 0.634 and R = 4.48 ■ 10" 4 . These values were 
derived from the numerical data using the first two moments 
of the distribution and Eqs. JT3}, (ITI)) . 



ready for f3Vo = 10 the whole distribution is well approx- 
imated by the stretched exponential with 7 « 0.90, Fig. 
IH1 Clearly, for ever higher barriers the transition kinet- 
ics becomes gradually single-exponential. This happens 
when the barrier height exceeds some characteristic value 
V c (a,T) which depends on a and temperature, cf. Fig. 
ITO1 Since there is no a precise threshold, the definition 
of V c is rather ambiguous. A working criterion for defin- 
ing V c can be, e.g., that the rate description is achieved 
within some error bound, e.g. 1% for deviation of 7 from 
unity. 

For f3V ^> (3V C , the overall escape kinetics is well- 
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FIG. 9: (Color online) Survival probability in one potential 
well versus time (in units of Tr b ') for (3Vo = 10. Using a 
fit to the survival probability, the whole distribution is well 
described by a stretched exponential with 7 ~ 0.90. The 
maximum likelihood fit to the numerical data yields but a 
slightly different value 7 ~ 0.85, cf. in Fig. 1101 indicating that 
the whole distribution is not properly stretched exponential. 
It is rather described by Eq. (|12[) with d = 0.335 and R = 
8.26 • 10" 6 . These values were derived from the numerical 
data using the first two moments of the distribution and Eqs. 

Ha), m. 



described by the non-Markovian rate theory. For a = 0.5 
it is very difficult to obtain a good statistics of transi- 
tions to find precisely this borderline. For the maximal 
value (3Vq — 12 used in our simulations the maximum 
likelihood fit with a stretched exponential (Weibull) dis- 
tribution yields 7 = 0.952 ±0.034. This value of f3V = 12 
provides an estimate for f3V c from below for a = 0.5. It 
approximately delimits the borderline between the appli- 
cability of NMRT for a = 0.5 and our treatment beyond 
it. The lower is a, the higher is borderline (3V c (a), and 
vice versa. For example, for a — 0.75, 0V C reduces to 
about (3V C ss 9. For these parameters, the maximum like- 
lihood fit of the numerical data with the single exponen- 
tial distribution yields the rate R = (2.864± 0.065) ■ 10" 5 
which only slightly deviates (about 3% of error) from the 
corresponding non-Markovian rate theory result R(n) ~ 
2.775 • IO" 5 . And for (3V = 10, the maximum likeli- 
hood fit yields R = (1.056 ± 0.034) • IO" 5 (a = 0.75) 
which almost agree within statistical errors with the non- 
Markovian rate theory result R(fi) = 1.021 ■ 10~ 5 . This 
provides a spectacular confirmation of both the non- 
Markovian rate theory for very high potential barriers, 
and the reliability of our numerics, as well as the physi- 
cal picture of anomalous escape developed in this work. 
On the contrary, for a < 0.5, e.g. a — 0.25, (3V C can be so 
high that the non-Markovian rate theory limit will never 
be reached for realistic barrier heights. Both our theoret- 
ical argumentation and the numerical results show that 
this borderline is fuzzy and and the rate description is 
restored gradually. The tendency in Fig. [TO] is, however, 
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FIG. 10: (Color online) Power exponent of a (single) stretched 
exponential fit to the overall residence time distribution ver- 
sus the barrier height for two different subdiffusion exponents 
a. The maximum likelihood approach is used to derive 7 from 
data. Notice that this 7 does not correspond to the asymp- 
totic 7 in the previous figures and text, but rather presents 
an average of ^(t) changing slowly in time. This subtle dif- 
ference diminishes when 7 approaches one. For a — 0.75, 
Co. 75 = 1.885; other parameters in the kernel approximation 
are the same as for a = 0.5. For Vb exceeding some borderline 
value V c (a, T), 7 tends to one and kinetics becomes gradually 
single exponential. 



obvious. 

The quasi-static disorder approximation cannot de- 
scribe quantitatively the numerical results for a broad 
range of parameters (rate disorder is yet dynamical, in 
spite of a quasi-infinite autocorrelation time [59j]). Nev- 
ertheless, it captures the essential physics (cf. Figs. 
I5I8|) and becomes ever better with increasing the bar- 
rier height, cf. Fig. [9J The agreement in this figure 
proves that our theory is essentially correct predicting 
the correct trend with increasing the barrier height, at 
least for a > 0.5. Indeed, with increasing the barrier 
height, or lowering the temperature the averaged escape 
time increases exponentially with /3Vb and, therefore, 
ever more slow noise components Q contributes to the 
non-Markovian rate and ever less such components con- 
tributes to fluctuation of this rate. For this reason, the 
root mean-squared amplitude of the slow (in our termi- 
nology) stochastic force £, s (t) gradually diminishes. For 
some characteristic [3V C , which clearly depends on a, it 
becomes negligible and the single-exponential kinetics is 
then approximately restored. The corresponding rate is 
given by the non-Markovian rate theory. 

The physical picture developed in this work is very 
different from the previous attempts in Refs. 0, to 
solve the problem of anomalous escape utilizing different 
approximations. The Ref. [H| focuses on the subdiffu- 
sive transmission through the parabolic barrier. It pre- 
dicts that asymptotic rate = — lim^oo dhi P(t)/dt 
always exists and is given precisely by the non-Markovian 
rate theory result in Eqs. (T5|)- (fTTT|) . This is clearly not cor- 



rect for V < V c (a,T). Strictly speaking, = always, 
even for V > V c (a,T) for a strictly power law memory 
kernel. However, for V 3> V c (a,T), the shape factor of 
Weibull distribution 7 equals approximately one, 7 sa I, 
and the rate description provides a good approximation. 
The higher is Vb, the better is this approximation, cf. 
Fig. M The Ref. focuses on the escape of a mass- 
less particle (m — > 0) from a parabolic potential well 
with a sharp cusp- like cutoff, utilizing the non-Markovian 
Fokker-Planck equation (NMFPE) of Refs. \M, 0- 
This NMFPE is exact for the parabolic potential, being 
but only approximate for a parabolic potential with cut- 
off. The better is the Gaussian approximation, the better 
should be the corresponding description, which implies 
high potential barriers (3Vq ^> 1. The theory in 0] can- 
not be compared directly with the present one (different 
potentials, zero mass particle in [4I(, expansion of the 
power law kernel into a finite sum of exponentials here), 
and the extrapolation of some main results in Ref. [4lj on 
a more realistic case here, would lead to the conclusions 
which are at odds with the present theory. In particular, 
the theory in 0] predicts (for a strict power law kernel, 
without inertial effects) that the escape kinetics is asymp- 
totically a power law, being only initially stretched expo- 
nential, and that the corresponding effective power law 
exponent tends exponentially to zero with increasing (3Vq. 
This means that the particle becomes strongly localized 
with increasing the barrier height, and the corresponding 
kinetics becomes ever more abnormal. On the contrary, 
the present theory predicts that the escape kinetics tends 
to a normal one, even if it decelerates dramatically. For 
a memory kernel with cutoff, the theory in Ref. [4 II ] pre- 
dicts that with increasing the barrier height the kinetics 
does become normal, when the memory cutoff becomes 
shorter than the mean escape time. This prediction con- 
cords with the present theory. The difference is however 
that the physical picture developed in this work suggests 
that the escape kinetics can be also approximately expo- 
nential when the memory cutoff largely exceeds the mean 
escape time. To conclude, the theory in this work is more 
physical. It overcomes the previous attempts to solve the 
very nontrivial problem of subdiffusive escape by taking 
a quite different road of multidimensional Markovian em- 
bedding and it is confirmed by numerics. 



The fact that the escape kinetics tends to a single- 
exponential with increasing the barrier height does not 
mean, however, that the diffusion becomes normal in the 
periodic potentials, as one might naively think in analogy 
with the CTRW theory. As a matter of fact, asymptot- 
ically such a diffusion cannot be faster that one in the 
absence of potential, i.e. (8x 2 (t)) oc t a . Therefore, we 
expect here new surprises. 



III. SUBDIFFUSION IN PERIODIC 
POTENTIALS 



We consider a common type washboard potential 
V(x) = — Vq cos(2itx/L) with the spatial period L. To 
study the influence of periodic potential on free subdif- 
fusion, it is convenient to scale now the time in the units 
of Tp — { y fiL 2 r\ a ) l l a , as in Fig. [21 which does not de- 
pend on the barrier height 2Vb- It takes time about Tp 
to subdiffuse freely over the distance about L. Indeed, 
the numerical simulations for /3Vo = 2 delivers a surprise 
indicating, see in Fig. [TTJ that the presence of peri- 
odic potential does not influence subdiffusion asymptoti- 
cally. This seems to agree with a theory in Refs. [39|, [42[ 
which, however, cannot be invoked directly because of it 
relates to a fully quantum case, where the tunneling ef- 
fects generally contribute. From this agreement we can, 
however, conclude that this surprising effect is certainly 
not of the quantum nature in the quantum case, but 
reflects the anti-persistent character of our viscoelastic 
subdiffusion which is purely classical. Namely, it is not 
diverging MRT, but extremally long-lived displacement 
(and velocity) anti-correlations which are responsible for 
the observed anomalous diffusion behavior in viscoelas- 
tic media. It must be emphasized, however, this asymp- 
totical regime is achieved through very long transients 
with a time-dependent a c s(t) gradually approaching a, 
cf. Fig. [11] This feature is to tally beyond any asymptot- 
ical analysis like one in Ref. [42| . The potential barrier 
height does generally matter and it strongly influences 
the whole time-course of diffusion. After a short ballistic 
stage followed by decaying coherent oscillations due to a 
combination of inertial and cage effects [61( in a poten- 
tial, the diffusion can look initially close to normal (as 
for f3Vo = 4 in Fig. ITT]) . This is due to a finite mean res- 
idence time in a potential well. However, it slows down 
and turns over into subdiffusion. The borderline of free 
subdiffusion cannot be crossed, cf. Fig. [TT] 

Very important is that the free viscoelastic subdiffu- 
sion is ergodic, in agreement with [64| . The results of 
the ensemble averaging with n = 10 4 particles coincide 
with the time averaging for a single particle done in ac- 
cordance with Eq. (|A5[) . see in Fig. [2a). However, in 
the periodic potential a strong deviation from the er- 
godic behavior takes place on the averaged time scale 
of the escape to the first neighboring potential wells, 
cf. Fig. QTJ This reflects anomalous escape kinetics 
as discussed above. Nevertheless, on a larger time-scale 
the single-trajectory averaging and the ensemble averag- 
ing yield again the same results. All this is in striking 
contrast with the CTRW-based subdiffusion, both free 
[lS EE IHHi and in periodic potentials [IS El]. In this 
respect, the benefit which subdiffusional search can bring 
for functioning the biological cells machinery [S H^] will 
not be questioned by the weak ergodicity breaking, as it 
might be in the case of CTRW-based subdiffusion. In- 
deed, the weak ergodicity breaking is related to a sponta- 
neous localization of CTRW-subdiffusing particles - i.e., 
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FIG. 11: (Color online) Diffusion in washboard potentials. 
The distance is measured in units of period L and the time in 
units of Tp = (/3L 2 r; a ) 1 ' /o! . a = 0.5 and the damping parame- 
ter rp = 773/Tbai = 10 with Tbai = L/vt and vt = ^Jk B T/m. 
Particles are initially localized at x = with Maxwellian dis- 
tributed velocities. The number of particles n = 10 4 is used in 
simulations (stochastic Heun algorithm, time step At = 10~ 4 ) 
for the ensemble averaging. Initially, the diffusional broaden- 
ing is always ballistic due to inertia effects. In potentials, 
this regime is followed by transient rattling oscillations due 
to a combination of the cage effect and the influence of the 
potential force. On the space scale 8x > L, diffusion be- 
comes ultraslow (5x 2 (t)) oc t a e£f(*) w jth a slowly changing 
a e ff(t) which depends on the potential height. Asymptoti- 
cally, a e ff(t — > 00) — ► 0.5. This universal asymptotics is 
almost reached for /3Vq = 2. 



a portion of them does not move at all (individual diffu- 
sion coefficient is close to zero), while other diffuse with 
an inhomogeneously distributed normal diffusion coeffi- 
cient [65], [66[ which depends on the total observation time 
T, even if the particles are totally identical. Numerically, 
or in a real experimental setup the ergodic behavior is 
achieved in the case of viscoelastic subdiffusion for very 
long T only, see also in [64j. So, to check the ergodicity 
for times until t = 10 3 we run a single trajectory for over- 
all T = 10 7 (this corresponds roughly to sampling over 
n = 10 4 copies, in analogy with the ensemble averaging). 
For a much smaller T, the difference between the time 
and ensemble averagings becomes sizeable. Hence, exper- 
imentally one can yet observe broadly distributed subd- 
iffusion coefficients, especially if both the particles and 
their environments are subjected to statistical variations 
However, differently from the CTRW subdiffusion, 
a particle will never get spontaneously trapped for the 
time of observation. This provides a true benchmark to 
distinguish between these two very different subdiffusion 
mechanisms experimentally. 



IV. SUMMARY 

Our main results are of profound importance for the 
anomalous diffusion and rate theory settling a long- 
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standing and controversial issue with conflicting results 
of different approaches, and different approximations. In 
particular, we prove that subdiffusion does not require 
principally a divergent mean residence time in a finite 
spatial domain, which makes it less anomalous when the 
anti-persistent, viscoelastic mechanism is at work. More- 
over, we substantiate the validity of the celebrated non- 
Markovian rate theory result (|5j)- (fTT)|) for very high po- 
tential barriers (j3V > 12 for a = 0.5 and (3V > (3V C « 9 
for a = 0.75 ) even for a strict power law memory kernel, 
where it was not expected to work because of an ultraslow 
relaxation within a potential well. However, for small 
a < 0.5 the corresponding borderline value (3V c {a) can 
be so high that this regime becomes practically unreach- 
able, at least for numerical simulations. Surprisingly, the 
non-Markovian theory result remains useful also for in- 
termediate barriers, 2ksT < Vq < V c (a,T), where it 
predicts the most probable logarithm of dwelling times. 
Here, the physics is well described by slowly fluctuating 
non-Markovian rates. For small barriers, (3Vq < 1, and 
for other models, e.g. when the bottom of potential well 
becomes more extended and flat, like the potential box 
in [&9( | , the fluctuating rate approach also loses its heuris- 
tic power. Then, the sluggish approach from the bottom 
to the barrier crossing region determines the transition 
kinetics. Even in the case of normal diffusion, different 
power-law kinetic regimes emerge 69] and the anomalous 
intrawell diffusion can change the corresponding power- 
law exponents, as modeled within the CTRW approach 
inRef. [H|. 

One of generic results is that the CTRW subdiffusion 
and the GLE subdiffusion are profoundly different, in 
spite of some superficial similarities. Subdiffusion in peri- 
odic potentials highlights the differences especially clear. 
Surprising is the finding that asymptotically the GLE 
subdiffusion is not sensitive to the barrier height, even 
if imposing a periodic potential does strongly affect the 
overall time-course of diffusion, and for a high potential 
barrier subdiffusion can look normal on a pretty long time 
interval. However, it slows down and asymptotically ap- 
proaches the borderline of free subdiffusion. Such subdif- 
fusion operates within a quite different (as compared with 
CTRW) physical mechanism based on the anti-persistent 
long-range correlations and not on the residence time dis- 
tributions with divergent mean. 

We believe that our results require to look anew on 
the theoretical interpretation of experimental subdiffu- 
sion results in biological applications, where the issue 
of ergodicity can be crucial. They provide some addi- 
tional theoretical support for the viscoelastic subdiffu- 
sion mechanism. A further detailed study is, however, 
necessary. To conclude, our work consolidates viscoelas- 
tic subdiffusion and fractional Brownian motion with the 
non-Markovian rate theory and fluctuating rate (dynam- 
ical disorder) approaches. It also agrees with the al- 
ready textbook view [see, e.g. in [7(| (pp. 380-382)], of 
the unusual kinetics as one with quasi-frozen and quasi- 
continuous conformational substates, as it was pioneered 



in biophysical applications by Austin el al.pl 



Acknowledgments 

Support of this work by the DFG-SFB-486 and by 
Volkswagen- Foundation (Germany), as well as useful dis- 
cussions with P. Talkner and P. Hanggi are gratefully 
acknowledged. 



APPENDIX A: EXACT SOLUTION OF THE 
POTENTIAL-FREE PROBLEM VERSUS MONTE 
CARLO SIMULATIONS 

In this Appendix, we discuss numerical errors by com- 
parison of the approximate results with the exact solution 
of subdiffusion problem in the absence of any potential. 
This exact solution is well-known [4p| . |47| . Assuming ini- 
tial velocities to be thermally distributed, it reads 



where 



(Sx 2 (t)) = 2v\ / H(t')dt', 
Jo 



H{t) = f K v (r)dT 
Jo 



(Al) 



(A2) 



is the integral of normalized equilibrium velocity auto- 
correlation function K v (t) := (v(t + r)v(t)) / (v 2 (0)) with 
(v 2 (0)) = v\ = ksT/ra. It has the Laplace-transform 



K v (s) 



1 



fj(s)/r 



(A3) 



Accordingly, the Laplace-transform of the coordinate 
variance, (Sx 2 (s)} := J °° exp(—st)(5x 2 (t))dt, is 



(5x 2 (s)) = 



2v\ 



s 2 [s + fj(s)/m] 



(A4) 



For the strict power-law kernel, Eq. (|A4p becomes 
(Sx 2 (s) = 2/[s 3 /r 2 + s 1+a ] with the distance measured 
in some arbitrary units L and time in the units of 
Tff = (L 2 r] a /kBT) 1 / a . This result can be inverted to the 
time domain in terms of the generalized Mittag-Lcfflcr 
functions, see e.g. in [6l|. However, both for the ex- 
act memory kernel and for its approximation in Eq. ([5]), 
it is convenient to invert Eq. (|A4[) numerically using the 
Gaver-Stehfest method [7l| with arbitrary numerical pre- 
cision, as done e.g. in Ref. [ll] for a different problem. 
The results thus obtained are numerically exact and the 
algorithm is very fast. They are compared against the 
results of the Monte-Carlo simulation of Eqs. J!])-® in 
Fig. [H In these simulations, we used both the ensemble 
averaging over n = 10 4 trajectories and a time averaging 
for single trajectory defined by 



(Sx 2 (t)) T = 



1 



T-t 



T-t 



[x(t + t')-x(t')] 2 dt' (A5) 
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for a very large time window T ^> t. The relative er- 
ror in Fig. [^b) is calculated as S — \(Sx^ xact (t)) — 
(fea P prox( i )>l/(^ X act(i)) I wh ere (Sx 2 applox (t)) is either 
the result of numerical solution of stochastic differen- 
tial equations (Monte Carlo, with 10 4 trajectories - noisy 
looking data), or the result of 16-exponential approxima- 
tion in Eq. (|A4[) . The agreement is confirming both for 
the used approximation of the memory kernel and for 
the quality of our stochastic simulations. The error in- 
troduced by the kernel approximation is mostly less than 
1%. The well-known phenomenon of logarithmic oscil- 
lations [l7| occurs within this error margin, and, there- 
fore, practically does not influence our stochastic numer- 
ics, which have a typical error of less than 3% (maximal 
5%). Notice that some damped oscillations in FigfTTlare 
of the inertial origin and have nothing in common with 



the logarithmic oscillations seen in Fig. [2jb). Moreover, 
for the octave scaling, b — 2, and for N — 64 such loga- 
rithmic oscillations are even not present in the variance 
behavior (not shown). The error margin in the variance 
behavior did not become, however, appreciably narrower. 
Therefore, such logarithmical oscillations practically do 
not matter for our numerics and the used N = 16 em- 
bedding computationally is even preferred. 

It worth to mention that the final time t = 10 in 
our simulations with n — 10 4 trajectories corresponds to 
about one week of computational time. Therefore, on our 
computers the theoretical limit of free normal diffusion 
for t > 10 12 for the used N = 16 embedding cannot 
be reached in principle. The presented data prove that 
our Markovian embedding is indeed both of a very good 
quality, and of practical use. 
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